Introduction

orchaRd 2.0 allows users to create pretty orchard plots that contain both confidence intervals (CIs) and prediction intervals (PIs) around point estimates for different categorical moderators, plots the effect size distribution over top such estimates and weights effect sizes by their precision (1/standard error, SE) or sample size. orchaRd takes a metafor object of class rma.mv or rma (Viechtbauer, 2010) and plots the results for the meta-analytic or meta-regression model.

orchaRd 2.0 now allows users to take meta-regression models with a multitude of moderator variables that are categorical or continuous along with heterogeneous residual variance models. This is achieved by allowing the user to produce marginalised mean estimates using the emmeans package. For plotting, orchaRd uses ggplot2 (Wickham, 2009), and as such, layers can be added directly to make plots customizable to the users needs. Also, the function orchard_plot uses the principles of bee swarm plots to visualise individual effect sizes (Eklund, 2012; see also Clarke & Sherrill-Mix, 2017; Sherrill-Mix & Clarke, 2017).

In addition to orchard plots, we have another plotting function caterpillars to draw caterpillar(s) plots. Furthermore, the package orchaRd comes with two functions that are useful for multilevel meta-analysis and meta-regression (i.e., rma.mv): 1) i2_ml to calculate the ‘total’ I2 along with I2 for each level (sensu Nakagawa S & Santos, 2012; see also A. M. Senior et al., 2016) and 2) r2_ml to calculate marginal and conditional R2 for mixed models (sensu S. Nakagawa & Schielzeth, 2013).

Citing orchaRd

To cite orchaRd in publications one can use the following reference:

Nakagawa, S. et al. 2021. The Orchard Plot: Cultivating the Forest Plot for Use in Ecology, Evolution and Beyond. Research Synthesis Methods, 12, 4–12.

Installation

To install orchaRd use the following code in R:

# install.packages('pacman')
rm(list = ls())
devtools::install_github("daniel1noble/orchaRd", force = TRUE)
pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, orchaRd, emmeans,
    ape, phytools, flextable)

Installation will make the primary functions accessible to users along with their help files. You will also need the tidyverse and metafor packages.

Examples of orchard and caterpillars plot

In this vignette we’ll walk the reader through a number of case studies and show you how you can create beautiful looking orchard plots. We overview three different case studies that make use of different effect sizes and moderators. The data sets associated with each case study come as part of the orchaRd package.

Example 1: Dietary Restriction and Lifespan

English & Uller (2016) performed a systematic review and meta-analysis on the effects of early life dietary restriction (a reduction in a major component of the diet without malnutrition; e.g. caloric restriction) on average age at death, using the standardised mean difference (often called d). They found that, across the whole data set, there was little evidence for an effect of dietary restriction on mean age at death. Here we’ll use the data set to first calculate the effect size measures and then fit an intercept only, meta-analytic model.

data(english)

# We need to calculate the effect sizes, in this case d
english <- escalc(measure = "SMD", n1i = NStartControl, sd1i = SD_C, m1i = MeanC,
    n2i = NStartExpt, sd2i = SD_E, m2i = MeanE, var.names = c("SMD", "vSMD"), data = english)

english_MA <- rma.mv(yi = SMD, V = vSMD, random = list(~1 | StudyNo, ~1 | EffectID),
    data = english)
summary(english_MA)
## 
## Multivariate Meta-Analysis Model (k = 77; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc  ​ 
## -44.5943   89.1887   95.1887  102.1809   95.5220   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0614  0.2478     21     no   StudyNo 
## sigma^2.2  0.0760  0.2756     77     no  EffectID 
## 
## Test for Heterogeneity:
## Q(df = 76) = 297.4722, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval    ci.lb   ci.ub   ​ 
##   0.0572  0.0729  0.7845  0.4327  -0.0856  0.2000    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We have fit a meta-analytic model and thus the only estimate we see is the overall effect size on the effects of caloric restriction on mean death across all studies examined. Now that we have fit our meta-analytic model we can get the confidence intervals and prediction intervals with a few functions in the orchaRd package. If one is interested in getting the table of results we can use the mod_results function. This will allow users to make nice tables of the results if needed. We can do that as follows:

model_results <- mod_results(english_MA, mod = "1", at = NULL, data = english, group = "StudyNo")
model_results
##      name estimate lowerCL upperCL lowerPR upperPR
## 1 Intrcpt    0.057  -0.086     0.2   -0.68     0.8

If we instead want to create an orchard plot and visualise the results we can do so quite simply as:

orchaRd::orchard_plot(english_MA, mod = "1", data = english, group = "StudyNo", xlab = "Standardised mean difference",
    transfm = "none")
Orchard plot of the impact caloric restriction using standardised mean difference

Figure 1. Orchard plot of the impact caloric restriction using standardised mean difference

In Figure 1, we simply add in the metafor model and it will create a default orchard plot. Alternatively, we could also add in the table of results. Also. some people may wish to present the heterogeneity index, I2 and, at the same time, change the colours of the ‘trunk and fruit’ by adding a few arguments.

# use i2_sn function to obtain the total I^2

I2 <- orchaRd::i2_ml(english_MA)

orchaRd::orchard_plot(model_results, mod="1", xlab = "Standardised mean difference") + 
  annotate(geom="text", x= 0.80, y= 1, label= paste0("italic(I)^{2} == ", round(I2[1],4)), # I want to use % but cannot do should come back
              color="black", parse = TRUE, size = 5) +
  scale_fill_manual(values="grey") +
  scale_colour_manual(values="grey")
Orchard plot of the impact caloric restriction using standardised mean difference but instead of using the metafor model, using the model results table

Figure 2. Orchard plot of the impact caloric restriction using standardised mean difference but instead of using the metafor model, using the model results table

Figure 2 and Figure 1 above show that overall estimate from a random-effects meta-analysis of 77 effect sizes is centered on zero, with a 95% CI that spans the line of no-effect. The prediction intervals clearly demonstrate the high level of heterogeneity (over 75%), with effects size less than -0.5 and greater than 0.5 predicted to be observed.

We can now draw a comparable caterpillar plot by using the function caterpillars.

# a caterpillar plot (not a caterpillars plot)
orchaRd::caterpillars(model_results, mod = "Int", xlab = "Standardised mean difference")
Catterpilar plot of the impact caloric restriction using standardised mean difference

Figure 3. Catterpilar plot of the impact caloric restriction using standardised mean difference

It is now interesting to compare how the same data set can be visualised different by look at Figure 2 and Figure 3. It is good that the caterpillar plot can show CIs for each effect size. An equivalent function was done with the size of fruit in the orchard plot.

In a subsequent publication, Alistair M. Senior, Nakagawa, Raubenheimer, Simpson, & Noble (2017) analysed this data set for effects of dietary-restriction on among-individual variation in the age at death using the log coefficient of variation ratio. A major prediction in both English & Uller (2016) and Alistair M. Senior et al. (2017) was that the type of manipulation, whether the study manipulated quality of food versus the quantity of food would be important. As such, we can fit a meta-regression model to test whether the moderator “Manipulation Type” ManipType impacts our results on the mean and variance

# First we need to calculate the lnCVR
english <- escalc(measure = "CVR", n1i = NStartControl, sd1i = SD_C, m1i = MeanC,
    n2i = NStartExpt, sd2i = SD_E, m2i = MeanE, var.names = c("lnCVR", "vlnCVR"),
    data = english)

# Now we can fit the meta-regression model (contrast)
english_MR0 <- rma.mv(yi = SMD, V = vSMD, mods = ~ManipType, random = list(~1 | StudyNo,
    ~1 | EffectID), data = english)
summary(english_MR0)
## 
## Multivariate Meta-Analysis Model (k = 77; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc  ​ 
## -44.2108   88.4216   96.4216  105.6916   96.9931   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0655  0.2560     21     no   StudyNo 
## sigma^2.2  0.0761  0.2758     77     no  EffectID 
## 
## Test for Residual Heterogeneity:
## QE(df = 75) = 295.5324, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.0414, p-val = 0.8388
## 
## Model Results:
## 
##                    estimate      se    zval    pval    ci.lb   ci.ub   ​ 
## intrcpt              0.0469  0.0924  0.5079  0.6115  -0.1342  0.2281    
## ManipTypeQuantity    0.0283  0.1390  0.2035  0.8388  -0.2442  0.3008    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Again, we can create a table of results
res2 <- orchaRd::mod_results(english_MR0, mod = "ManipType", data = english, group = "StudyNo")
res2
##       name estimate lowerCL upperCL lowerPR upperPR
## 1  Quality    0.047   -0.13    0.23   -0.71    0.81
## 2 Quantity    0.075   -0.14    0.30   -0.69    0.84
# we fit a meta-regression (contrast) for the log coefficient of variation
senior_MR0 <- rma.mv(yi = lnCVR, V = vlnCVR, mods = ~ManipType, random = list(~1 |
    StudyNo, ~1 | EffectID), data = english)
summary(senior_MR0)
## 
## Multivariate Meta-Analysis Model (k = 77; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc  ​ 
## -32.0740   64.1481   72.1481   81.4180   72.7195   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0275  0.1657     21     no   StudyNo 
## sigma^2.2  0.0470  0.2169     77     no  EffectID 
## 
## Test for Residual Heterogeneity:
## QE(df = 75) = 215.7242, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 1.3308, p-val = 0.2487
## 
## Model Results:
## 
##                    estimate      se     zval    pval    ci.lb   ci.ub   ​ 
## intrcpt             -0.1310  0.0678  -1.9333  0.0532  -0.2639  0.0018  . 
## ManipTypeQuantity    0.1217  0.1055   1.1536  0.2487  -0.0851  0.3285    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# creating a table of results
res3 <- orchaRd::mod_results(senior_MR0, mod = "ManipType", data = english, group = "StudyNo",
    at = list(ManipType = "Quality"), subset = TRUE)
res3
##      name estimate lowerCL upperCL lowerPR upperPR
## 1 Quality    -0.13   -0.26  0.0018   -0.68    0.42
# We can now plot SMD and lnCVR beside each other and compare the results
p1 <- orchaRd::orchard_plot(res2, mod = "ManipType", data = english, group = "StudyNo",
    xlab = "Standardised mean difference")

p2 <- orchaRd::orchard_plot(res3, mod = "ManipType", data = english, group = "StudyNo",
    xlab = "log(CV ratio) (lnCVR)")

p1/p2
Orchard plot of diet qualities impact on SMD (top) and log coefficient of variation (bottom)

Figure 4. Orchard plot of diet qualities impact on SMD (top) and log coefficient of variation (bottom)

Our orchard plot for the log coefficient of variation demonstrates that, while restrictions on dietary quality and quantity do not affect the average age at death (top of Figure 4), among-individual variation may be altered by quality restrictions (bottom of Figure 4). The effect is negative suggesting that the coefficient of variation in the control group is lower than that in the treatment group, and the 95% CI does not span zero. Again though, the effect is heterogeneous; a substantial number of positive effects are still predicted.

Let’s draw comparable caterpillars plots to compare these two results again.

# We can now plot SMD and lnCVR beside each other and compare the results
p1c <- orchaRd::caterpillars(res2, mod = "ManipType", data = english, group = "StudyNo",
    xlab = "Standardised mean difference")

p2c <- orchaRd::caterpillars(res3, mod = "ManipType", data = english, group = "StudyNo",
    xlab = "log(CV ratio) (lnCVR)")

p1c/p2c
Catterpilar plot of diet qualities impact on the log coefficient of variation (lnCVR)

Figure 5. Catterpilar plot of diet qualities impact on the log coefficient of variation (lnCVR)

Now compare between Figure 4 and Figure 5. We think they both look nice.

Example 2: Predation and Invertebrate Community

Eklof et al. (2012) evaluated the effects of predation on benthic invertebrate communities. Using the log response ratio they quantified differences in abundance and/or biomass of gastropods and Amphipods in groups with and without predation in an experimental setting.

Here again, we can create orchard plots of the model results, but we’ll show how a few simple things can be modified. Again, we can fit the meta-analytic model first:

data(eklof)

# Calculate the effect size
eklof <- escalc(measure = "ROM", n1i = N_control, sd1i = SD_control, m1i = mean_control,
    n2i = N_treatment, sd2i = SD_treatment, m2i = mean_treatment, var.names = c("lnRR",
        "vlnRR"), data = eklof)

# Add the observation level factor
eklof$Datapoint <- as.factor(seq(1, dim(eklof)[1], 1))

# Also, we can get the sample size, which we can use for weighting if we would
# like
eklof$N <- rowSums(eklof[, c("N_control", "N_treatment")])

# fit a meta-regression with the intercept (contrast)
eklof_MR0 <- rma.mv(yi = lnRR, V = vlnRR, mods = ~Grazer.type, random = list(~1 |
    ExptID, ~1 | Datapoint), data = eklof)

summary(eklof_MR0)
## 
## Multivariate Meta-Analysis Model (k = 34; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc  ​ 
## -51.9349  103.8698  111.8698  117.7328  113.3513   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed     factor 
## sigma^2.1  0.5075  0.7124     18     no     ExptID 
## sigma^2.2  0.6703  0.8187     34     no  Datapoint 
## 
## Test for Residual Heterogeneity:
## QE(df = 32) = 195.9943, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.8891, p-val = 0.3457
## 
## Model Results:
## 
##                       estimate      se     zval    pval    ci.lb    ci.ub    ​ 
## intrcpt                -0.8095  0.3099  -2.6119  0.0090  -1.4170  -0.2021  ** 
## Grazer.typegastropod    0.3285  0.3484   0.9429  0.3457  -0.3544   1.0114     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Above we have fit a meta-regression model using “Grazer Type” as a moderator which is predicted to explain variation in log response ratios. We can demonstrate a few simple changes users can make, but we note here that users can make far more complex changes down the line if needed, but we’ll save explaining these more complex changes for the last example. The first is the angle at which the y-axis labels are positioned (bottom of Figure 6):

f <- orchaRd::mod_results(eklof_MR0, mod = "Grazer.type", data = eklof, group = "ExptID")

p3 <- orchaRd::orchard_plot(eklof_MR0, mod = "Grazer.type", data = eklof, group = "ExptID",
    xlab = "log(Response ratio) (lnRR)", transfm = "none")

p4 <- orchaRd::orchard_plot(eklof_MR0, mod = "Grazer.type", data = eklof, group = "ExptID",
    xlab = "log(Response ratio) (lnRR)", transfm = "none", angle = 45, g = FALSE)

p3/p4
Orchard plots of the effects of predation on benthic invertebrate communities compared using the log response ratio. Top panel is the default plot and bottom panel a plot containing changes in label axes

Figure 6. Orchard plots of the effects of predation on benthic invertebrate communities compared using the log response ratio. Top panel is the default plot and bottom panel a plot containing changes in label axes

The other thing we can change is the type of scaling we wish to use. Let’s say we are interested in scaling the effect size by the total sample size of the study we use a vector of N, sample size (bottom of Figure 7:

p5 <- orchaRd::orchard_plot(eklof_MR0, mod = "Grazer.type", data = eklof, group = "ExptID",
    xlab = "log(Response ratio) (lnRR)")

p6 <- orchaRd::orchard_plot(eklof_MR0, mod = "Grazer.type", data = eklof, group = "ExptID",
    xlab = "log(Response ratio) (lnRR)", angle = 45, N = "N", g = FALSE)

p5/p6
Orchard plots of the effects of predation on benthic invertebrate communities compared using the log response ratio. Top panel is the default plot and bottom panel a plot containing changes in label axes and scaling with sample size instead of precision

Figure 7. Orchard plots of the effects of predation on benthic invertebrate communities compared using the log response ratio. Top panel is the default plot and bottom panel a plot containing changes in label axes and scaling with sample size instead of precision

Overall, our orchard plot shows the results of a re-analysis of their data. The estimated mean effects are negative for both gastropods and amphipods suggesting that mean abundance/biomass in the control group is lower than in the treatment groups, although the effect is largest, and is statistically significant, for amphipods. In both cases the prediction intervals reveal the extent of heterogeneity, with positive effects predicted to be observed.

Again, we can also draw the caterpillars-plot couterpart for this data set (bottom of Figure 8):

eklof_MR_result <- orchaRd::mod_results(eklof_MR0, mod = "Grazer.type", data = eklof,
    group = "ExptID")
orchaRd::caterpillars(eklof_MR_result, mod = "Grazer.type", data = eklof, group = "ExptID",
    xlab = "log(Response ratio) (lnRR)")
Caterpillars plot of the effects of predation on benthic invertebrate communities compared using the log response ratio

Figure 8. Caterpillars plot of the effects of predation on benthic invertebrate communities compared using the log response ratio

Like last time, we compare between Figure 7 and Figure 8. Again, they both look quite informative.

Example 3: Maternal-Offspring Morphological Correlations

Finally, we also look at the case discussed by Lim, Senior, & Nakagawa (2014), who meta-analysed the strength of correlation between maternal and offspring size within-species, across a very wide range of taxa. They found, that typically, there is a moderate positive correlation between maternal size and offspring size within species (i.e. larger mothers have larger offspring). However, they also found evidence for relatively strong phylogenetic effects suggesting that the strength of the association is dependent on evolutionary lineage.

Here we have used an orchard plot to represent the results obtained when meta-analysing the data from Lim et al. (2014) by taxonomic Phylum.

data(lim)

# Add in the sampling variance
lim$vi <- (1/sqrt(lim$N - 3))^2

# Lets fit a meta-regression - I will do Article non-independence. The
# phylogenetic model found phylogenetic effects, however, instead we could fit
# Phylum as a fixed effect and explore them with an Orchard Plot
lim_MR <- metafor::rma.mv(yi = yi, V = vi, mods = ~Phylum - 1, random = list(~1 |
    Article, ~1 | Datapoint), data = lim)
summary(lim_MR)
## 
## Multivariate Meta-Analysis Model (k = 357; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc  ​ 
## -97.6524  195.3049  213.3049  248.0263  213.8343   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed     factor 
## sigma^2.1  0.0411  0.2029    220     no    Article 
## sigma^2.2  0.0309  0.1757    357     no  Datapoint 
## 
## Test for Residual Heterogeneity:
## QE(df = 350) = 1912.9637, p-val < .0001
## 
## Test of Moderators (coefficients 1:7):
## QM(df = 7) = 356.6775, p-val < .0001
## 
## Model Results:
## 
##                        estimate      se     zval    pval    ci.lb   ci.ub     ​ 
## PhylumArthropoda         0.2690  0.0509   5.2829  <.0001   0.1692  0.3687  *** 
## PhylumChordata           0.3908  0.0224  17.4190  <.0001   0.3468  0.4347  *** 
## PhylumEchinodermata      0.8582  0.3902   2.1992  0.0279   0.0934  1.6230    * 
## PhylumMollusca           0.4867  0.1275   3.8175  0.0001   0.2368  0.7366  *** 
## PhylumNematoda           0.4477  0.3054   1.4658  0.1427  -0.1509  1.0463      
## PhylumPlatyhelminthes    0.4935  0.2745   1.7980  0.0722  -0.0444  1.0314    . 
## PhylumRotifera           0.4722  0.3021   1.5634  0.1180  -0.1198  1.0642      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Noe we can plot a default orchard plot, scaling each effect size by N. Also, because we are using Zr, we can use transfm = “tanh” and it will do the conversions for us:

# Plot the meta-regression model
orchaRd::orchard_plot(lim_MR, mod = "Phylum", group = "Article", data = lim, xlab = "Correlation coefficient",
    alpha = 0.5, transfm = "tanh", angle = 45, N = "N", cb = FALSE)
Orchard plot of the the strength of correlation between maternal and offspring size within-species

Figure 9. Orchard plot of the the strength of correlation between maternal and offspring size within-species

Now that we have Figure 9, we can do some small changes to make it pretty and add some more details, such as between study heterogeneity and R2 for phylum (we use ). Currently, the cb argument is “FALSE”, we can change this to “TRUE” to use colour blind friendly colours. Additionally, because we are using ggplot we can add element to the figure to make it look pretty.

# Lets add the R2 on the figure as well as little modifications:
R2 <- orchaRd::r2_ml(lim_MR)

# We can draw an orchard plot with R2
orchaRd::orchard_plot(lim_MR, mod = "Phylum", group = "Article", data = lim, xlab = "Correlation coefficient (r)",
    alpha = 0.5, transfm = "tanh", angle = 45, N = "N", cb = TRUE) + theme(legend.position = c(0.05,
    0.99), legend.justification = c(0, 1), legend.key.size = unit(1, "mm")) + theme(legend.direction = "horizontal",
    legend.title = element_text(size = 8), legend.text = element_text(size = 10)) +
    annotate(geom = "text", x = 0.8, y = 0.7, label = paste0("italic(R)^{2} == ",
        round(R2[1], 4) * 100), color = "black", parse = TRUE, size = 4)
Orchard plot of the the strength of correlation between maternal and offspring size within-species

Figure 10. Orchard plot of the the strength of correlation between maternal and offspring size within-species

As in Figure 10, new elements can be added to the orchard_plot to modify it as one sees fit. It will overwrite existing elements. From our orchard plots above, it is clear that the analysis is dominated by data from Chordates and Arthropods, with the other Phyla being much more poorly represented. Second, there is a difference between the strength of a typical correlation within these two well represented groups (the correlation is stronger in Chordates), which arguably would explain the phylogenetic signals detected by Lim et al. (2014). Lastly, although there are differences within the typical correlation between Chordates and Arthropods, there remains a large overlap in predicted range of individual effect sizes; individual species within Phyla are still highly variable.

Finally, we make a comparable caterpillars plot.

# Plot the meta-regression model
lim_MR_results <- orchaRd::mod_results(lim_MR, mod = "Phylum", group = "Article",
    data = lim)
orchaRd::caterpillars(lim_MR_results, mod = "Phylum", group = "Article", data = lim,
    xlab = "Correlation coefficient", transfm = "tanh", g = FALSE)
Caterpillars plot of the the strength of correlation between maternal and offspring size within-species

Figure 11. Caterpillars plot of the the strength of correlation between maternal and offspring size within-species

Compared to the orchard plot (Figure 10, Figure 11) does not look as elegant, unlike the previous comparisons. This is because the groups with small sample sizes are packed at the top of the figure.

From Figure 10 we can see that many taxonimic groups have very little data. We may want instead to just focus the readers attention to means of groups with sufficient data. The new mod_results function allows us to do this using the maginalised means approach. We can do that as follows:

# Plot the meta-regression model
lim_MR_results <- orchaRd::mod_results(lim_MR, mod = "Phylum", group = "Article",
    data = lim, at = list(Phylum = c("Chordata", "Arthropoda", "Mollusca")), subset = TRUE)
orchaRd::orchard_plot(lim_MR_results, data = lim, xlab = "Correlation coefficient",
    transfm = "tanh", g = TRUE, angle = 45)
Caterpillars plot of the the strength of correlation between maternal and offspring size within-species

Figure 12. Caterpillars plot of the the strength of correlation between maternal and offspring size within-species

We can now see that Figure 12 limits what it plots to only the three levels relevant for presentation.

Orchard Plots with Meta-regression Models: Marginalised means and heteroscedastic residual variances

In orchaRd 2.0 we’ve made some major updates to the types of models that can be used with the package. Meta-analyses in ecology and evolution often fit a multitude of moderators and random effects within a single model. This is even desirable in many cases to deal with high heterogeneity and disentangle the independent contributions of moderators to changes in mean effect size. The challenge with the previous version of the package was that it was difficult to accommodate these models. However, the emmeans package is an excellent solution to the problem. It will produce marginsalised mean estimates for a given level of moderator – producing mean estimates that average over all other moderators in the model. This includes models with categorical and continuous moderators.

An additional issue that is of importance with meta-regression models is the need to deal with heteroscedastic variance across levels of a moderator. Homogeneity of variance assumptions are often violated in meta-analysis. orchaRd 2.0 allows the user to fit models that explicitly estimate different residual variances reducing the probability of type I errors.

In this section, we show users: 1) how to fit these models using metafor and 2) how these models interface with the orchaRd package to produce pretty plots that capture the model predictions.

Example 4: Meta-regression and heterogeneous variance models

Meta-regression models with single categorical moderators are extremely common. Meta-analyst’s want to estimate the mean effect size within levels of a moderator and test whether these means are significantly different from each other. One inherent assumption to fitting these models is that we are assuming that the within study (or residual) variance within these groups is the same (i.e., “homogeneity of variance assumption”). This assumption will likely be violated in many situations in meta-regression models. Unless one conducts a sub-group analysis for each level of the moderator (although this is less powerful and makes it hard to test whether levels are significantly different) different residual variances need to be explicitly modeled in metafor. Fortunately, this can be done. orchaRd 2.0 now allows plotting of these models with their corresponding CI’s and PI’s.

We’ll work with a meta-analysis by O’Dea, Lagisz, Hendry, & Nakagawa (2019). In this meta-analysis, the authors were interested in testing whether developmental temperature in fish had an effect on mean and variance in a host of phenotypic traits. They used experimental studies manipulating developmental temperatures using the log response ratio (lnrr) and the log coefficient of variation (lncvr) as their effect sizes. We’ll only focus on lnrr here for demonstration purposes. Importantly, they categorised effect sizes as having come from one of 4 trait categories: Behaviour, Life-history, Morphology and Physiology. They were interested in whether developmental temperature impact trait types differently. We’ll first explore how to fit a meta-regression model that explicitly estimates residual variance in each of these trait levels.

First, load the data.

# Load the dataset that comes with orchaRd
data(fish)

# Subset data for demonstration purposes.
warm_dat <- fish

To fit a heterogeneous variance model in metafor we need to use some additional arguments and specify the random effect structure differently (i.e., specify inner argument). To do this, we need to add the struc = HCS to the argument list, which specifies the variance structure of an ~inner|outer formula of a particular random effect. In this case, we’re assuming a heteroscedastic compound symmetry structure. To avoid estimating the correlation / covariance between levels we need to also add the rho argument and set this to 0. This will assume that the covariance between the levels is set to 0 (although it could be estimated).

Once we have these additional arguments we need to modify our random effect structure. Note that metafor does not estimate a residual / within study variance. We need to explicitly add an observation level random effect, in this case es_ID. If we have more than one effect size / study (which is common) then we want to have this random effect in our models.

model_het <- metafor::rma.mv(yi = lnrr, V = lnrr_vi, mods = ~trait.type, method = "REML",
    test = "t", random = list(~1 | group_ID, ~1 + trait.type | es_ID), rho = 0, struc = "HCS",
    data = warm_dat, control = list(optimizer = "optim", optmethod = "Nelder-Mead"))
model_het
## 
## Multivariate Meta-Analysis Model (k = 410; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2    0.0143  0.1196     84     no  group_ID 
## 
## outer factor: es_ID      (nlvls = 410)
## inner factor: trait.type (nlvls = 4)
## 
##             estim    sqrt  k.lvl  fixed         level 
## tau^2.1    0.2144  0.4630      4     no     behaviour 
## tau^2.2    0.2121  0.4605      4     no  life-history 
## tau^2.3    0.0232  0.1522    323     no    morphology 
## tau^2.4    0.0320  0.1790     79     no    physiology 
## rho        0.0000                   yes               
## 
## Test for Residual Heterogeneity:
## QE(df = 406) = 39294.5728, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 406) = 3.0890, p-val = 0.0271
## 
## Model Results:
## 
##                         estimate      se     tval   df    pval    ci.lb   ci.ub​ 
## intrcpt                   0.1848  0.2713   0.6811  406  0.4962  -0.3485  0.7180 
## trait.typelife-history    0.5576  0.3701   1.5068  406  0.1326  -0.1699  1.2851 
## trait.typemorphology     -0.1846  0.2714  -0.6801  406  0.4968  -0.7181  0.3489 
## trait.typephysiology     -0.1712  0.2730  -0.6272  406  0.5309  -0.7079  0.3655 
##  
## intrcpt 
## trait.typelife-history 
## trait.typemorphology 
## trait.typephysiology 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can see from above that we now have four variances; one estimated for each level of the trait.type moderator. Now that we have fit our model and estimated the relevant parameters we can produce the orchard plot using the marginal_means function:

het_model <- orchaRd::mod_results(model_het, group = "group_ID", mod = "trait.type",
    data = warm_dat)

orchaRd::orchard_plot(het_model, xlab = "lnRR", data = warm_dat)
Orchard plot with heteroscedastic error for trait typess

Figure 13. Orchard plot with heteroscedastic error for trait typess

We can now see from Figure 13 that the variance in effects for the morphological and physiological trait categories are smaller than for life-history and behavioural trait categories; likely because of the very small sample sizes in the latter trait categories.

As a second example, we can turn to Pottier, Burke, Drobniak, Lagisz, & Nakagawa (2021) who did a meta-analysis comparing acclimation differences between males and females. Pottier et al. (2021) found that aquatic and terrestrial systems differed in their variability of effects. Here, we’ll also show how phylogenetic correlation matrices can be incorporated within the models and work with orchaRd.

# Clear working space
rm(list = ls())

# Load the data
data("pottier")

# Load the phylogenetic correlation matrix
data("phylo_matrix")

# Load the sampling variance matrix
data("VCV_dARR")

Now that we have all the data and matrices loaded, we can fit the multi-level mete-regression model.

##### Heteroscedasticity modeled at the effect size level
mod.habitat_het <- rma.mv(yi = dARR, V = VCV_dARR, mods = ~habitat, method = "REML",
    test = "t", dfs = "contain", random = list(~1 | species_ID, ~1 | phylogeny, ~habitat |
        es_ID), struct = "HCS", rho = 0, R = list(phylogeny = phylo_matrix), data = pottier)
mod.habitat_het
## 
## Multivariate Meta-Analysis Model (k = 1089; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed      factor    R 
## sigma^2.1  0.0090  0.0948    138     no  species_ID   no 
## sigma^2.2  0.0129  0.1136    138     no   phylogeny  yes 
## 
## outer factor: es_ID   (nlvls = 1089)
## inner factor: habitat (nlvls = 2)
## 
##             estim    sqrt  k.lvl  fixed        level 
## tau^2.1    0.0716  0.2677    929     no      aquatic 
## tau^2.2    0.0082  0.0907    160     no  terrestrial 
## rho        0.0000                   yes              
## 
## Test for Residual Heterogeneity:
## QE(df = 1087) = 56246.7341, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 136) = 9.2383, p-val = 0.0028
## 
## Model Results:
## 
##                     estimate      se     tval   df    pval    ci.lb    ci.ub​ 
## intrcpt               0.2087  0.0655   3.1856  136  0.0018   0.0792   0.3383 
## habitatterrestrial   -0.1492  0.0491  -3.0394  136  0.0028  -0.2463  -0.0521 
##  
## intrcpt             ** 
## habitatterrestrial  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now we can produce an orchaRd plot using this model to show the unequal variance in the two habitat types.

orchard_plot(mod.habitat_het, data = pottier, group = "species_ID", mod = "habitat",
    xlab = "dARR", angle = 45)
Orchard plot with heteroscedastic error for habitat type

Figure 14. Orchard plot with heteroscedastic error for habitat type

flextable::flextable(mod_results(mod.habitat_het, data = pottier, group = "species_ID",
    mod = "habitat")$mod_table)

We might also think about heterogeneous variance existing at multiple random effect levels. These models require a lot of data, and need to be carefully thought about, but we orchard 2.0 can also handle these models.

##### Heteroscedasticity modeled at the effect size level
mod.habitat_het2 <- rma.mv(yi = dARR, V = VCV_dARR, mods = ~habitat, method = "REML",
    test = "t", dfs = "contain", random = list(~habitat | species_ID, ~1 | phylogeny,
        ~habitat | es_ID), struct = "HCS", rho = 0, phi = 0, R = list(phylogeny = phylo_matrix),
    data = pottier)
mod.habitat_het2
## 
## Multivariate Meta-Analysis Model (k = 1089; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed     factor    R 
## sigma^2    0.0098  0.0989    138     no  phylogeny  yes 
## 
## outer factor: species_ID (nlvls = 138)
## inner factor: habitat    (nlvls = 2)
## 
##             estim    sqrt  k.lvl  fixed        level 
## tau^2.1    0.0116  0.1077    929     no      aquatic 
## tau^2.2    0.0022  0.0469    160     no  terrestrial 
## rho        0.0000                   yes              
## 
## outer factor: es_ID   (nlvls = 1089)
## inner factor: habitat (nlvls = 2)
## 
##               estim    sqrt  k.lvl  fixed        level 
## gamma^2.1    0.0715  0.2673    929     no      aquatic 
## gamma^2.2    0.0084  0.0918    160     no  terrestrial 
## phi          0.0000                   yes              
## 
## Test for Residual Heterogeneity:
## QE(df = 1087) = 56246.7341, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 136) = 14.6883, p-val = 0.0002
## 
## Model Results:
## 
##                     estimate      se     tval   df    pval    ci.lb    ci.ub​ 
## intrcpt               0.2098  0.0581   3.6096  136  0.0004   0.0949   0.3247 
## habitatterrestrial   -0.1606  0.0419  -3.8325  136  0.0002  -0.2434  -0.0777 
##  
## intrcpt             *** 
## habitatterrestrial  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
orchard_plot(mod.habitat_het2, data = pottier, group = "species_ID", mod = "habitat",
    xlab = "dARR", angle = 45)
Orchard plot with heteroscedastic error for habitat type at the residual and species level

Figure 15. Orchard plot with heteroscedastic error for habitat type at the residual and species level

flextable::flextable(mod_results(mod.habitat_het2, data = pottier, group = "species_ID",
    mod = "habitat")$mod_table)

Example 5: Meta-regression and marginalised mean effect sizes

We’ll now demonstrate how multi-moderator meta-regression models can be used with orchaRd 2.0, by plotting marginalised mean effects within categories of a single moderator. Marginalised means are also extremely useful to improve interpretation of overall model results.

To demonstrate how to do this we’ll again turn to O’Dea et al. (2019) as our example data set. We’ll first load the data set, if you haven’t already done so:

# Load the dataset that comes with orchaRd
data(fish)

# Subset data for demonstration purposes.
warm_dat <- fish

Now that we have the data set assume we are interested in modelling multiple moderators that we hypothesize will explain effect size variance. We want to provide a model that will allow us to understand the conditional (independent) effects of each moderator on mean effect size. We know that experimental_design and the temperature difference between treatments (deg_diff) are expected to impact mean effect size. We also expect experiments that last longer to have stronger effects (treat_end_days). We’ll now extend our model from Section 4 above to include these moderators so we can look at the effects of each controlling for each other.

# Fit the metaregerssion model
model <- metafor::rma.mv(yi = lnrr, V = lnrr_vi, mods = ~experimental_design + trait.type +
    deg_dif + treat_end_days, method = "REML", test = "t", random = list(~1 | group_ID,
    ~1 + trait.type | es_ID), rho = 0, struc = "HCS", data = warm_dat, control = list(optimizer = "optim",
    optmethod = "Nelder-Mead"))

model
## 
## Multivariate Meta-Analysis Model (k = 410; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2    0.0143  0.1196     84     no  group_ID 
## 
## outer factor: es_ID      (nlvls = 410)
## inner factor: trait.type (nlvls = 4)
## 
##             estim    sqrt  k.lvl  fixed         level 
## tau^2.1    0.2145  0.4631      4     no     behaviour 
## tau^2.2    0.2147  0.4634      4     no  life-history 
## tau^2.3    0.0226  0.1503    323     no    morphology 
## tau^2.4    0.0320  0.1790     79     no    physiology 
## rho        0.0000                   yes               
## 
## Test for Residual Heterogeneity:
## QE(df = 402) = 38537.0410, p-val < .0001
## 
## Test of Moderators (coefficients 2:8):
## F(df1 = 7, df2 = 402) = 2.4871, p-val = 0.0165
## 
## Model Results:
## 
##                                           estimate      se     tval   df​ 
## intrcpt                                     0.1637  0.2726   0.6003  402 
## experimental_designsplit pooled families    0.0339  0.0503   0.6728  402 
## experimental_designsplit single family      0.0654  0.0687   0.9513  402 
## trait.typelife-history                      0.5532  0.3721   1.4870  402 
## trait.typemorphology                       -0.1847  0.2726  -0.6775  402 
## trait.typephysiology                       -0.1689  0.2745  -0.6153  402 
## deg_dif                                    -0.0091  0.0051  -1.7891  402 
## treat_end_days                              0.0007  0.0003   2.2183  402 
##                                             pval    ci.lb   ci.ub 
## intrcpt                                   0.5486  -0.3723  0.6997    
## experimental_designsplit pooled families  0.5015  -0.0651  0.1328    
## experimental_designsplit single family    0.3420  -0.0697  0.2005    
## trait.typelife-history                    0.1378  -0.1782  1.2847    
## trait.typemorphology                      0.4985  -0.7206  0.3512    
## trait.typephysiology                      0.5387  -0.7085  0.3707    
## deg_dif                                   0.0744  -0.0190  0.0009  . 
## treat_end_days                            0.0271   0.0001  0.0014  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Our model is complicated. It not only estimates the effects of trait type, but type of experimental design, degree difference and treatment duration. But, effects of moderators are now conditional on each other. What happens if we just want to plot mean effects for trait type? We can do this using marginalised means. We can see that there’s not much going on for experimental design, so we will marginalize the means across each level of this moderator. Degree difference and treatment end days are continuous, however. To deal with continuous moderators we need to specify the levels we want predictions to be made. We can do that using the marginal_means function:

HetModel <- orchaRd::mod_results(model, group = "group_ID", mod = "trait.type", at = list(deg_dif = c(5,
    10, 15)), by = "deg_dif", weights = "prop", data = warm_dat)

orchaRd::orchard_plot(HetModel, xlab = "lnRR", data = warm_dat, angle = 45, g = FALSE,
    legend.pos = "top.left") + theme(legend.direction = "vertical")
Orchard plot of marginalised mean effect sizes for 5, 10, and 15 degree Celcius temperture differences

Figure 16. Orchard plot of marginalised mean effect sizes for 5, 10, and 15 degree Celcius temperture differences

We can see the orchard plot changes. Here we have a condition label. This specifies the degree difference for mean predicted. Note that, these means effectively average over all levels of experimental design and all possible values of treat_end_days. We have also made some small changes to the plot. First, we set g = FALSE which will suppress the total groups on the plot (e.g. studies), moved the position of legends and angled the labels.

We can also specify the specific value of treat_end_day if the user would like to be more precise. We can do that as follows:

HetModel2 <- orchaRd::mod_results(model, group = "group_ID", mod = "trait.type",
    at = list(deg_dif = c(5, 10, 15), treat_end_days = c(10)), by = "deg_dif", weights = "prop",
    data = warm_dat)

orchaRd::orchard_plot(HetModel2, xlab = "lnRR", data = warm_dat, angle = 45, g = FALSE,
    legend.pos = "top.left") + theme(legend.direction = "vertical")
Orchard plot of marginalised mean effect sizes for 5, 10, and 15 degree Celcius temperture differences

Figure 17. Orchard plot of marginalised mean effect sizes for 5, 10, and 15 degree Celcius temperture differences

Now the orchard plot depicts the mean effect sizes (plus 95% CI and PI) for each trait category, at 5, 10 and 15 \(^{\circ}\)C when the treatment duration is 10 days.

References

Clarke, E., & Sherrill-Mix, S. (2017). Ggbeeswarm: Categorical scatter (violin point) plots. R Package Version 0.6.0.
Eklof, J. S., Alsterberg, C., Havenhand, J. N., Sundback, K., Wood, H. L., & Gamfeldt, L. (2012). Experimental climate change weakens the insurance effect of biodiversity. Ecology Letters, 15, 864–872.
Eklund, A. (2012). Beeswarm: The bee swarm plot, an alternative to stripchart. R Package Version 0.1.
English, S., & Uller, T. (2016). Does early-life diet affect longevity? A meta-analysis across experimental studies. Biology Letters, 12, 20160291.
Lim, J. N., Senior, A. M., & Nakagawa, S. (2014). Heterogeneity in individual quality and reproductive trade-offs within species. Evolution, 68(8), 2306–18. doi: 10.1111/evo.12446.
Nakagawa, S, & Santos, E. S. A. (2012). Methodological issues and advances in biological meta-analysis. Evolutionary Ecology, 26(5), 1253–1274.
Nakagawa, S., & Schielzeth, H. (2013). A general and simple method for obtaining R2 from generalized linear mixed‐effects models. Methods in Ecology and Evolution, 4(2), 133–142.
O’Dea, R. E., Lagisz, M., Hendry, A. P., & Nakagawa, S. (2019). Developmental temperature affects phenotypic means and variability: A meta-analysis of fish data. Fish and Fisheries, 20, 1005–1022.
Pottier, P., Burke, S., Drobniak, S. M., Lagisz, M., & Nakagawa, S. (2021). Sexual (in) equality? A meta‐analysis of sex differences in thermal acclimation capacity across ectotherms. Functional Ecology, in press, https://doi.org/10.1111/1365–2435.13899.
Senior, A. M., Grueber, C. E., Kamiya, T., Lagisz, M., O’Dwyer, K., Santos, E. S. A., & S, N. (2016). Heterogeneity in ecological and evolutionary meta-analyses: Its magnitudes and implications. Ecology, 97(12), 3293–3299.
Senior, Alistair M., Nakagawa, S., Raubenheimer, D., Simpson, S. J., & Noble, D. W. A. (2017). Dietary restriction increases variability in longevity. Biology Letters, 13(3), 20170057.
Sherrill-Mix, S., & Clarke, E. (2017). Plot categorical data using quasirandom noise and density estimates. R Package Version 0.4.5.
Viechtbauer, W. (2010). Conducting meta-analyses in r with the metafor package. Journal of Statistical Software, 363, 1–48.
Wickham, H. (2009). ggplot2: Elegant graphics for data analysis. New York ; London: Springer.